#include "defaults.h"
#include "nurbs.h"
#include "utils.h"
#include "serialize.h"
#include "primitivemanager.h"
#include "triangle/triangle.h"
#include <sstream>
#include <boost/timer.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>

using namespace std;

extern PrimitiveManager* primitiveManager;

std::string Filename(const std::string& str) {
	using namespace boost::posix_time;
	using namespace boost::gregorian;
	ptime now = second_clock::local_time();
	date today = now.date();
	std::stringstream ss1(stringstream::out);
	ss1 << str << to_iso_string(now) << ".";
	return ss1.str();
}

std::string GeneratePath(const std::string& strPrefix) {
	using namespace boost::posix_time;
	using namespace boost::gregorian;
	ptime now = second_clock::local_time();
	date today = now.date();
	std::stringstream ss1(stringstream::out);
	ss1 << strPrefix << to_iso_string(now);
	return ss1.str();
}

std::string FileExtension(const std::string& filename) {
	std::string str = filename;
	const size_t extensionLength = str.length() - str.rfind('.') - 1;

	str = str.substr(str.size() - extensionLength);

	return str;
}

vector<nfloat> DetermineParameters(const Knotvector& U, size_t n, size_t pu, nfloat uStart, nfloat uEnd, size_t samplesPatch) {
	vector<nfloat> U2;
	U2.push_back(uStart);
	
	size_t uspanStart = (size_t)math::spline::FindSpan(n, pu, uStart, U);
	size_t uspanEnd = (size_t)math::spline::FindSpan(n, pu, uEnd, U);

	for(size_t i = uspanStart; i < uspanEnd; i++)
		U2.push_back(U[i + 1]);
	U2.push_back(uEnd);
	
	vector<nfloat> Unew;
	for(size_t i = 0; i < U2.size() - 1; i++) {
		nfloat a = U2[i];
		nfloat b = U2[i + 1];
		if(std::abs(b - a) > 0.001) {
			for(size_t j = 0; j < samplesPatch; j++) {
				nfloat t = a + ((nfloat)j / samplesPatch) * (b - a);
				Unew.push_back(t);
			}
		}
	}
	Unew.push_back(uEnd);

	return Unew;
}

////////////////////////////////////////////////////////////////////////////////////////////////
size_t GetIndex(const Vec4n& P, map<string, size_t>& mapPosition, vector<Vec4n>& position, const nfloat eps = 0.00000001) {
	stringstream ss;

	ss << fixed;
 	Vec3n P2 = (1.0 / P[3]) * math::V3n(P[0], P[1], P[2]);
	for(size_t i = 0; i < 3; i++) if(std::abs(P2[i]) < eps) P2[i] = 0.0;
	ss << P2[0] << " " << P2[1] << " " << P2[2];

	if(mapPosition[ss.str()] == 0) {
		size_t i = position.size() + 1;
		mapPosition[ss.str()] = i;
		position.push_back(P);
		return i - 1;
	}
	
	return mapPosition[ss.str()] - 1;	
}

template<size_t D>
size_t GetIndex(const math::Vec<nfloat, D>& P, std::vector< math::Vec<nfloat, D> >& position, const nfloat eps = 0.00000001) {
	for(size_t i = 0; i < position.size(); i++) {
		if(length(position[i] - P) < eps)
			return i;
	}
	position.push_back(P);
	return position.size() - 1;
}

void ExtractPatches(const Surface& surface, map<string, size_t>& mapPosition, vector<Vec4n>& position, 
					vector< boost::multi_array<size_t, 2> >& element, vector< pair<Vec2n, Vec2n> >& domain, vector< pair<Knotvector, Knotvector> >& knotvector) {
	const Knotvector& U = surface.GetKnotvectorU();
	const Knotvector& V = surface.GetKnotvectorV();
	size_t p = surface.GetDegreeU();
	size_t q = surface.GetDegreeV();
	
	for(size_t i = 0; i < surface.Size0() - p; i++) {
		for(size_t j = 0; j < surface.Size1() - q; j++) {	
			//domain
			Vec2n uDomain = math::V2n(U[i + p], U[i + p + 1]);
			Vec2n vDomain = math::V2n(V[j + q], V[j + q + 1]);
			const nfloat eps = 0.0000001;
			if(std::abs(uDomain[1] - uDomain[0]) > eps && std::abs(vDomain[1] - vDomain[0]) > eps) {
				domain.push_back(pair<Vec2n, Vec2n>(uDomain, vDomain));
				
				//coefficients				
				element.push_back(boost::multi_array<size_t, 2>(boost::extents[p + 1][q + 1]));
				boost::multi_array<size_t, 2>& c = element.back();
				for(size_t i0 = 0; i0 <= p; i0++) {
					for(size_t j0 = 0; j0 <= q; j0++) {
						const Vec4n& P = surface(i + i0, j + j0);
// 						c[i0][j0] = GetIndex(P, mapPosition, position);
						c[i0][j0] = GetIndex<4>(P, position);
					}
				}
				//knots
				Knotvector Ui, Vi;
				for(size_t k = 0; k <= 2 * p + 1; k++) Ui.push_back(U[i + k]);
				for(size_t k = 0; k <= 2 * q + 1; k++) Vi.push_back(V[j + k]);
				knotvector.push_back(pair<Knotvector, Knotvector>(Ui, Vi));
			}
		}
	}
}

NurbsSurfaceMesh* GenerateSurface(const std::vector<Surface>& surface, bool useWeights) {
	vector<Vec4n> position4;
	vector< boost::multi_array<size_t, 2> > element;
	vector< std::pair<Vec2n, Vec2n> > domain;
	vector< std::pair<Knotvector, Knotvector> > knotvector;
	std::vector<size_t> sid;

	size_t oldsize = 0;
	map<string, size_t> mapPosition;
 	for(size_t i = 0; i < surface.size(); i++) {
		ExtractPatches(surface[i], mapPosition, position4, element, domain, knotvector);
		for(size_t j = 0; j < element.size() - oldsize; j++)
			sid.push_back(i);
		oldsize = element.size();
	}

	vector<Vec3n> position(position4.size());
	vector<nfloat> weight(position.size());
	for(size_t i = 0; i < position.size(); i++) {
		const Vec4n& P = position4[i];
		weight[i] = P[3];
		position[i] = (1.0 / weight[i]) * math::V3n(P[0], P[1], P[2]);
	}
	PRINTVAR(position.size());

	std::vector< std::pair<std::string, Variant> > vertexdata, nurbsdata;
	vertexdata.push_back(std::pair<std::string, Variant>("position", position));
 	if(useWeights)
		vertexdata.push_back(std::pair<std::string, Variant>("weight", weight));
	nurbsdata.push_back(std::pair<std::string, Variant>("sid", sid));

	const size_t p = surface[0].GetDegreeU();
	const size_t q = surface[0].GetDegreeV();

	return new NurbsSurfaceMesh(position.size(), vertexdata, nurbsdata, element, knotvector, p, q);
}

NurbsSurfaceMesh* GenerateSurface(const Surface& surface, bool useWeights) {
	std::vector<Surface> s;
	s.push_back(surface);
	return GenerateSurface(s, useWeights);
}

Surface GenerateTorus() {
	const size_t p = 2;
	const size_t q = 2;
	const size_t n = 9;
	const size_t m = 9;
	nfloat tu[12] = {0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4};
	nfloat tv[12] = {0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4};

	nfloat c[9][9][3] = {{{5, 0,-1}, {6, 0,-1}, {6, 0, 0}, {6, 0, 1}, {5, 0, 1}, {4, 0, 1}, {4, 0, 0}, {4, 0,-1}, {5, 0,-1}},
						 {{5, 5,-1}, {6, 6,-1}, {6, 6, 0}, {6, 6, 1}, {5, 5, 1}, {4, 4, 1}, {4, 4, 0}, {4, 4,-1}, {5, 5,-1}},
						 {{0, 5,-1}, {0, 6,-1}, {0, 6, 0}, {0, 6, 1}, {0, 5, 1}, {0, 4, 1}, {0, 4, 0}, {0, 4,-1}, {0, 5,-1}},
						 {{-5, 5,-1}, {-6, 6,-1}, {-6, 6, 0}, {-6, 6, 1}, {-5, 5, 1}, {-4, 4, 1}, {-4, 4, 0}, {-4, 4,-1}, {-5, 5,-1}},
						 {{-5, 0,-1}, {-6, 0,-1}, {-6, 0, 0}, {-6, 0, 1}, {-5, 0, 1}, {-4, 0, 1}, {-4, 0, 0}, {-4, 0,-1}, {-5, 0,-1}},
						 {{-5,-5,-1}, {-6,-6,-1}, {-6,-6, 0}, {-6,-6, 1}, {-5,-5, 1}, {-4,-4, 1}, {-4,-4, 0}, {-4,-4,-1}, {-5,-5,-1}},
						 {{0,-5,-1}, {0,-6,-1}, {0,-6, 0}, {0,-6, 1}, {0,-5, 1}, {0,-4, 1}, {0,-4, 0}, {0,-4,-1}, {0,-5,-1}},
						 {{5,-5,-1}, {6,-6,-1}, {6,-6, 0}, {6,-6, 1}, {5,-5, 1}, {4,-4, 1}, {4,-4, 0}, {4,-4,-1}, {5,-5,-1}},
						 {{5, 0,-1}, {6, 0,-1}, {6, 0, 0}, {6, 0, 1}, {5, 0, 1}, {4, 0, 1}, {4, 0, 0}, {4, 0,-1}, {5, 0,-1}}};

	nfloat w[9][9] = {{1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1},
					  {0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548},
					  {1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1},
					  {0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548},
					  {1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1},
					  {0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548},
					  {1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1},
					  {0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548, 0.5, 0.707106781186548},
					  {1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1, 0.707106781186548, 1}};

	Knotvector U;
	for(size_t i = 0; i < 12; i++)
		U.push_back(0.25 * tu[i]);

	Knotvector V;
	for(size_t i = 0; i < 12; i++)
		V.push_back(0.25 * tv[i]);

	Array2V4n Pw(boost::extents[n][m]);
	for(size_t i = 0; i < n; i++) {
		for(size_t j = 0; j < m; j++) {
			Vec3n x = math::V3n(c[i][j][0], c[i][j][1], c[i][j][2]);

			Pw[i][j] = w[i][j] * math::V4n(x, 1.0);

		}
	}

	return Surface(Pw, p, q, U, V);
}

Surface Square(size_t N, size_t M, nfloat a, nfloat b, size_t p, size_t q, bool identity = false) {
	size_t n = N - 1; //number of control points per slice
	size_t m = M - 1; //number of slices
	vector< vector<Vec4n> > Pw(N);
	nfloat xoffset = 0.0;
	nfloat yoffset = 0.0;
	Knotvector U = math::spline::GenerateOpenKnotvector<nfloat>(p, N);
	Knotvector V = math::spline::GenerateOpenKnotvector<nfloat>(q, M);

	size_t multiplicity = 1;

	Knotvector U2 = U;
	for(size_t i = p + 1; i < U.size() - p - 1; i += multiplicity) {
		nfloat t = 0.0;
		for(size_t j = 0; j < multiplicity; j++)
			t += U[i + j];
		t /= (nfloat)multiplicity;

		for(size_t j = 0; j < multiplicity; j++)
			U2[i + j] = t;
	}
	U = U2;

	Knotvector V2 = V;
	for(size_t i = q + 1; i < V.size() - q - 1; i += multiplicity) {
		nfloat t = 0.0;
		for(size_t j = 0; j < multiplicity; j++)
			t += V[i + j];
		t /= (nfloat)multiplicity;

		for(size_t j = 0; j < multiplicity; j++)
			V2[i + j] = t;
	}
	V = V2;

	nfloat w = 1.0;
	if(!identity) {
		for(size_t i = 0; i < N; i++) {
			nfloat x = (nfloat)i / (N - 1);
			for(size_t j = 0; j < M; j++) {
				nfloat y = (nfloat)j / (M - 1);
				Vec4n p;
				p = w * math::V4n(a * x, b * y, 0.0, 1.0);
				Pw[i].push_back(p);
			}
		}
	} else {
		for(size_t i = 0; i < N; i++) {
			nfloat x = math::spline::Node<nfloat>(U, i, p);
			for(size_t j = 0; j < M; j++) {
				nfloat y = math::spline::Node<nfloat>(V, j, q);
				Vec4n P;
				P = w * math::V4n(a * x, b * y, 0.0, 1.0);
				Pw[i].push_back(P);
			}
		}
	}

	return Surface(Pw, p, q, U, V);
}

Surface DiskPolar(nfloat scale) {
	std::vector< std::vector<Vec4n> > cPw(2);
	for(int i = 0; i < 2; i++)
		cPw[i].resize(9);
	
	cPw[0][0] = math::V4n(1.0, 0.0, 0.0, 1.0);
	cPw[0][1] = math::V4n(std::sqrt(0.5), std::sqrt(0.5), 0.0, std::sqrt(0.5));
	cPw[0][2] = math::V4n(0.0, 1.0, 0.0, 1.0);
	cPw[0][3] = math::V4n(-std::sqrt(0.5), std::sqrt(0.5), 0.0, std::sqrt(0.5));
	cPw[0][4] = math::V4n(-1.0, 0.0, 0.0, 1.0);
	cPw[0][5] = math::V4n(-std::sqrt(0.5), -std::sqrt(0.5), 0.0, std::sqrt(0.5));
	cPw[0][6] = math::V4n(0.0, -1.0, 0.0, 1.0);
	cPw[0][7] = math::V4n(std::sqrt(0.5), -std::sqrt(0.5), 0.0, std::sqrt(0.5));
	cPw[0][8] = math::V4n(1.0, 0.0, 0.0, 1.0);

	cPw[1][0] = math::V4n(0.0, 0.0, 0.0, 1.0);
	cPw[1][1] = math::V4n(0.0, 0.0, 0.0, std::sqrt(0.5));
	cPw[1][2] = math::V4n(0.0, 0.0, 0.0, 1.0);
	cPw[1][3] = math::V4n(0.0, 0.0, 0.0, std::sqrt(0.5));
	cPw[1][4] = math::V4n(0.0, 0.0, 0.0, 1.0);
	cPw[1][5] = math::V4n(0.0, 0.0, 0.0, std::sqrt(0.5));
	cPw[1][6] = math::V4n(0.0, 0.0, 0.0, 1.0);
	cPw[1][7] = math::V4n(0.0, 0.0, 0.0, std::sqrt(0.5));
	cPw[1][8] = math::V4n(0.0, 0.0, 0.0, 1.0);
	Knotvector cU;
	cU.push_back(0.0);
	cU.push_back(0.0);
	cU.push_back(1.0);
	cU.push_back(1.0);
	Knotvector cV;
	cV.push_back(0.0);
	cV.push_back(0.0);
	cV.push_back(0.0);
	cV.push_back(0.25);
	cV.push_back(0.25);
	cV.push_back(0.5);
	cV.push_back(0.5);
	cV.push_back(0.75);
	cV.push_back(0.75);
	cV.push_back(1.0);
	cV.push_back(1.0);
	cV.push_back(1.0);
	Surface polar(cPw, 1, 2, cU, cV);

	if(std::abs(scale - 1.0) < 0.00001) {
		polar.DegreeElevateU(1);
		return polar;
	}

	Mat4x4n R;
	R.Identity();
	R.RotateZ(45.0);
	for(size_t i = 0; i < polar.Size0(); i++) {
		for(size_t j = 0; j < polar.Size1(); j++) {
			Vec4n& P = polar(i, j);
			P = R * P;


		}
	}

	vector<Knotvector> X(2);
	X[0].push_back(scale);
	polar.Refine(X);

	PRINTVAR(polar.Size0());
	PRINTVAR(polar.Size1());
	vector< vector<Vec4n> > Pw(polar.Size0() - 1);
	for(size_t i = 0; i < Pw.size(); i++) {
		Pw[i] = vector<Vec4n>(polar.Size1());
		for(size_t j = 0; j < Pw[i].size(); j++) {
			Pw[i][j] = polar(i, j);
		}
	}

	Knotvector U = math::spline::GenerateOpenKnotvector<nfloat>(1, 2);
	Knotvector V = polar.GetKnotvectorV();

	Surface polar2(Pw, 1, 2, U, V);
	polar2.DegreeElevateU(1);
	return polar2;
}

Curve Circle() {
	std::vector<Vec4n> Pw(9);
	Pw[0] = math::V4n(1.0, 0.0, 0.0, 1.0);
	Pw[1] = math::V4n(std::sqrt(0.5), std::sqrt(0.5), 0.0, std::sqrt(0.5));
	Pw[2] = math::V4n(0.0, 1.0, 0.0, 1.0);
	Pw[3] = math::V4n(-std::sqrt(0.5), std::sqrt(0.5), 0.0, std::sqrt(0.5));
	Pw[4] = math::V4n(-1.0, 0.0, 0.0, 1.0);
	Pw[5] = math::V4n(-std::sqrt(0.5), -std::sqrt(0.5), 0.0, std::sqrt(0.5));
	Pw[6] = math::V4n(0.0, -1.0, 0.0, 1.0);
	Pw[7] = math::V4n(std::sqrt(0.5), -std::sqrt(0.5), 0.0, std::sqrt(0.5));
	Pw[8] = math::V4n(1.0, 0.0, 0.0, 1.0);
	
	Knotvector U;
	U.push_back(0.0);
	U.push_back(0.0);
	U.push_back(0.0);
	U.push_back(0.25);
	U.push_back(0.25);
	U.push_back(0.5);
	U.push_back(0.5);
	U.push_back(0.75);
	U.push_back(0.75);
	U.push_back(1.0);
	U.push_back(1.0);
	U.push_back(1.0);
			
	return Curve(Pw, 2, U);
}

Surface PeriodicCylindricalSurface(const nfloat radius, const nfloat height, const size_t n, const size_t m, const size_t p, const size_t q) {
	const Knotvector U = math::spline::GenerateOpenKnotvector<nfloat>(p, n);
	const Knotvector V = math::spline::GeneratePeriodicKnotvector<nfloat>(q, m);

	Curve curve = Circle();

	std::vector<Vec3n> Q;
	for(size_t j = 0; j < m; j++) {
		const nfloat t = (nfloat)j / m;
		Q.push_back(radius * curve.Evaluate(t));
	}
	for(size_t j = 0; j < q; j++)
		Q.push_back(Q[j]);

	std::reverse(Q.begin(), Q.end());

	Array2V4n Pw(boost::extents[n][Q.size()]);
	for(size_t i = 0; i < n; i++) {
		const nfloat h = height * math::spline::Node<nfloat>(U, i, p);

		for(size_t j = 0; j < Q.size(); j++) {
			Pw[i][j] = math::V4n(Q[j][0], Q[j][1], h, 1.0);
		}
	}
	
	return Surface(Pw, p, q, U, V);
}
